% Gets value function by iteration on discrete grid
VDIFF=inf;                                                 % reset difference in value function
Viter=0;                                                   % reset counter

Vf = beta/(1-beta)*PAYOFFMAT;

while VDIFF>sqrt(eps) && Viter<20000                       % iterate to convergence of V 
  Viter=Viter+1;                                           % increment counter

% Time t+1 values
% [EV, pstar,exitflag] = get_pstar(Vf, Pgrid, mu);                 % quadratic value function interpolation 
% if exitflag==0
%     sound(sin(1:300));
%     zero = NaN(size(in)); % tell solver to check another point
%     disp('V function interpolation not working properly')
%     return
% end

logitprobMAT = logitprob(Vf, noise, wbar);
EV = sum(logitprobMAT.*Vf);
% relentropy = ones(nump,1)*(log(nump) + sum(logitprobMAT.*log(max(logitprobMAT,eps^.5))));
relentropy = [];
exitflag = 1;
P = logitprobMAT; 
pstar = Pgrid'*logitprobMAT;
% P = Pmatrix(pstar,Pgrid,pstep);  
% logitprobMAT = P;
% relentropy = [];

adjcost = cost(model, menucost, noise, relentropy);
D = ones(nump,1)*EV - R'*T2'*Vf - wbar*adjcost;                 % D is the expected gain from adjustment (can be negative)
Dopt = ones(nump,1)*max(Vf) - R'*T2'*Vf - wbar*adjcost;         % Dopt is the gain from adjustment to rational optimum
lambda = adjustment(model, D, alpha, alphapos, menucost, lbar,omega,Pgrid, pstar);  %lambda = adjustment(model, D/wbar, alpha, menucost, lbar);
    
CONTINVALUE = Vf + D.*lambda;   
  
% Time t values
iterV = PAYOFFMAT + beta*R'*T2'*CONTINVALUE;  % iterV is current payoff plus disc. continuation value
%iterV = PAYOFFMAT + beta*(R'*T2'*Vf+D.*lambda);  % iterV is current payoff plus disc. continuation value


VdifMin = min(min(iterV-Vf));                             % minimum difference 
VdifMax = max(max(iterV-Vf));                             % maximum difference
VDIFF = VdifMax-VdifMin;                                  % distance between max and min difference
Vf = iterV;                                               % updating V
end
Vf = Vf + beta/(1-beta)*(VdifMax+VdifMin)/2;                % vertical shift of value function once shape has converged

if VDIFF>sqrt(eps) || any(isnan(Vf(:)))
    warning('Vf has not converged!');
    sound(sin(1:300));
    zero = NaN; % tell solver to check another point
    return;
end